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Abstract. It has been suggested that the nonlinear Schrodinger-Newton equation 
might approximate the coupling of quantum mechanics with gravitation, particularly 
in the context of the M0ller-Rosenfeld semiclassical theory. Numerical results for the 
spherically symmetric, time-dependent, single-particle case are presented, clarifying 
and extending previous work on the subject. It is found that, for a particle mass 
greater than 1.14(?i^/(G'ct))^/'^, a wave packet of width a partially "collapses" to 
a groundstate solution found by Moroz, Penrose, and Tod, with excess probability 
dispersing away. However, for a mass less than 1.14(?i^/(G'(j))^/'^, the entire wave 
packet appears to spread like a free particle, albeit more slowly. It is argued that, on 
some scales (lower than the Planck scale), this theory predicts significant deviation 
from conventional (linear) quantum mechanics. However, owing to the difficulty of 
controlling quantum coherence on the one hand, and the weakness of gravity on the 
other, definitive experimental falsification poses a technologically formidable challenge. 



PACS numbers: 04.62. +v,03.65.-w,04.20.-q 
1. Introduction 

Noting tliere is no direct evidence that the gravitational field is quantized, Carlip has 
suggested the possibility that it is not [I]. The suggestion is not new, and Carlip points 
specifically to a theory put forth by Rosenfeld [2j and independently by M0ller [3], that 
the gravitational field retains its classical nature while coupling to an expectation value 
of quantized matter: 

Gafe = 87r(r,fe). (1) 

Although a superficially similar equation has been discussed in the context of a 
semiclassical approximation to quantum gravity [1] , M0ller and Rosenfeld proposed that 
Eq. ([T]) is exact. Despite its unconventional character, Carlip has argued that the theory 
has not yet been ruled out, and with Salzman has worked towards an experimentally 
testable prediction O [6l [1]. The intent of the present work is to continue this effort 
with more accurate numerical results. 
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In the nonrelativistic limit, {Tab) = {^\Tab\^) where \I/ is a state function that 
evolves according to Schrodinger's equation [UIH], 

= - E^V^v^(X,t) -Gf:m,$(x„t)vI/(X,t), (2) 

i=l * 1=1 

and the potential $ is given by the Newtonian limit of Eq. ([T]), assuming source 
particles, 

N 

V2$(a;,t) =47rG' / d3^X'|^(X', t)p ^ m,5(x - x;), (3) 

i=i 

and X = (xi, . . . , Xat). Eq. ([2]) coupled with Eq. ([3]) is known as the Schrodinger-Newton 
system, which can be written as a single equation by solving Eq. ([3]) with the Green's 
function: 

= -f:I^Vf*(X.<)-G f: ,„,,„,/ <;-.Y'^5;^*(X,()(4) 

Conventionally, a Schrodinger potential excludes nonlinear self-interaction. In 
violating this rule, Eq. (jl]) represents a fundamental alteration of quantum mechanics. 
First, because it includes self- interaction, the potential is present for any massive particle 
irrespective of the presence of other particles or any external potential. There are no 
free particles in this theory. Second, it is nonlinear even for a single particle. Thus, 
arbitrary linear superpositions are no longer permitted (although they will often serve 
as good approximations when the mass is small enough that the nonlinear term can be 
neglected) . 

It should be emphasized that this Schrodinger-Newton system differs markedly from 
the conventional coupling of the Schrodinger equation with a Newtonian or Coulomb- 
like potential, as noted in [9j. It is true that the Schrodinger-Newton equation shares 
terms in common with the standard mean-field approximation of such a potential. For 
example, substituting the Hartree ansatz ^ = 'i/'(xi,t) ■ • ■il){'Xj\r,t) into the conventional 
system gives: 

ih = - - — Vi?/'i(xi,t) -G^minij / d x- ^-7— ?/'i(xi, t) 

dt 2mi p{ J |xj-x^| 

+ Gm,^/d^a;: 'f^^;'^i|% .(x.,t) (5) 

J I Xj X j I 

(Usually the sum in Eq. (|5]) excludes the j = i term; here the j = i term is included 
in the sum while its negation is separately added, for emphasis.) Applying the same 
ansatz to the Schrodinger-Newton system yields: 

ih = -- — VjV^i(xi,t) -G^minij / d x - ^-^-^'4x^,^06) 

ot Irrii ~^ J |Xj — x^l 

These equations are similar but Eq. (jH]) includes self-interaction, whereas Eq. 
excludes it. This distinction is also true of the Hartree- Fock approach, with an ansatz 
consisting of the Slater determinant and a potential corrected for fermion statistics. In 
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this case the final term in Eq. is replaced by a sum called the "exchange term" 
which, again, cancels the self-interaction term in the conventional case. 

This distinguishing feature is most obvious for a single-particle system. In this 
case, the standard mean-field potential vanishes, leaving the equation for a free particle. 
Meanwhile the Schrodinger-Newton system retains the self-interaction term: 

^n^-^ = -^VM^,t)-Gm^ /rfVM^^(x,t). (7) 
at 2m J |x — x'l 

In the many-particle case, the significant physical difference between the 
Schrodinger-Newton and conventional Hartree or Hartree-Fock systems is made more 
evident by transforming to relative and center-of-mass coordinates, R = I^ili ^jXi/M, 
M = Y^f^iTTii, Tj = Xj — xat, j = 1,...,A^ — 1, and using the ansatz = 
'i/'_R(R, t)^/'ri (fi) ■ ■ ■ ^r]v_i (riv-i, (or the Slater determinant). In the conventional 
system, the potential can be expressed entirely in terms of relative coordinates. 
Separation of variables then leaves the center-of-mass equation identical in form to that 
of a free-particle with mass M and position R. In the Schrodinger-Newton system, 
on the other hand, the potential cannot be expressed entirely in terms of relative 
coordinates, because of the appearance of the kernel |xAr — x'^|~^. In this case the 
center-of-mass receives self-interaction, and separation of variables results in a center- 
of-mass equation identical in form to the nonlinear, single-particle Schrodinger-Newton 
equation (Eq. ([7j)), with mass M and position R. 

There have been various efforts to investigate physical consequences of this 
nonlinearity. Diosi argued that it would result in soliton-like behavior, which he 
predicted would effectively localize macroscopic objects [8J. He further estimated the 
width and energy of the single-particle groundstate. Independently Penrose conjectured 
that quantum gravity requires gravitationally induced collapse of the wave function, and 
suggested the Schrodinger-Newton equation might approximate this process [IDl [HI Il2] • 
This motivated a numerical study by Moroz, Penrose, and Tod (MPT), which yielded 
stationary solutions in the single particle case, including a groundstate consistent with 
Diosi's estimate [T2] . 

Subsequent numerical studies focused on the time- dependent, single-particle 
Schrodinger-Newton equation. The following were among the features observed of the 
solution set: 

• Some solutions tended to contract in radius, for the most part evolving towards a 
narrower wave function. This was observed by Harrison et al. [I3l [lH [15], Salzman 
and Carlip [5l[6l[T], and Giulini and Grossardt [16]. Guzman et al. also made similar 
observations of Schrodinger-Newton solutions albeit in a different physical context, 
that of hypothetical scalar fields in astrophysical systems [171 llHl [IS] . Harrison et 
al. identified the narrower wave function as a scaled-down MPT groundstate. 

• These same solutions fragmented slightly, releasing pulses carrying a small fraction 
of their total probability to infinity, while the rest of their amplitude remained more 
localized. This was observed by Harrison et al. and Guzman et al. 



Schrddinger-Newton "collapse" of the wave function 



4 



• Below a critical mass, solutions were found which seemed to expand, in their 
entirety, without bound, much like free particles but with slower velocities. This 
was observed by Salzman and Carlip, and Giulini and Grossardt, although the 
critical mass found by the latter was six orders of magnitude greater than that 
found by the former. 

The present study qualitatively demonstrates all of the above features. 
Quantitatively, the results presented here are consistent with Harrison et al. in regards 
to the final state towards which some solutions evolve, to within a few percent, and with 
Giulini and Grossardt in regards to the value of the critical mass below which solutions 
expand towards spatial infinity, to within ~ 20%. (Note the present work was originally 
submitted before the Giulini and Grossardt preprint appeared.) 

We remark in passing that the term "collapse" , which is used to describe aspects of 
Schrodinger-Newton evolution here and in [5l|6l[T], should be qualified. The term is used 
because in some respects the dynamics resemble the traditional conception of collapse 
of the wave function: for an amount of matter approaching macroscopic scales, the wave 
function tends to localize in position and approach an energy eigenstate. However, unlike 
some traditional notions about wave function collapse, Schrodinger-Newton evolution is 
continuous and unitary. 

The present investigation will focus on the numerical evolution of a single particle 
with the same initial data used by Salzman and Carlip: 

V>(r,0) = (7ra2)-3/V^V(2.^) (g) 

This paper is organized as follows: In Sec. 2 the numerical methods are explained. In 
Sec. 3 the main numerical results are given. And in Sec. 4 physical implications and 
the possibility for empirical verification are discussed. 



2. Methods 



As in [5] for a single-particle source and spherical symmetry the Poisson equation 
(Eq. ([3])) can be solved using the appropriate Green's function, 

1 



yielding 



47r max(r, r') ' 



G{f,r')\ij{7^,t)\^dV' 



(9) 



- AnGm 



-1 



r'^\tlj{r',t)\''dr' + I r'\^{r' ,t)\''dr' 



(10) 



In the present work this is integrated with a 4th order accurate Newton-Cotes algorithm 
(i.e. the error is (9(Ar^) where Ar is the interval between adjacent points of the 
discretized radius). Fig. [1] demonstrates its accuracy by comparing its result for the 
initial $ to the analytic solution, Gm^erf(r/cr)/r, obtained from Eq. (fTOj) with the 
initial ip of Eq. 1^. 
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Figure 1. Relative error for the initial computation of $ is shown, at two different 
resolutions. The scaling is such that for the expected 4th order accuracy the two curves 
should coincide. 
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As in [15], Schrodinger's equation can be simplified slightly in spherical symmetry 
by the change of variables u = rip: 

dtu = i — d'^u — i—^u (11) 
2m n 

Here the spatial derivative on the right hand side is computed with a 6th order accurate 
stencil. To provide a boundary condition appropriate for finite differencing near the 
origin, u{r) is assumed to be effectively an odd function (as is consistent with the initial 
data and Eq. ( Ilip ). Time- integration is then performed with a 4th order accurate 
Runge-Kutta algorithm (that is, 4th order in At). Note that stability requires that 
At oc Ar^ Thus the time-integration is 8th order in Ar, and the dominant error, outside 
of the Poisson-integration, should be 6th order in Ar. Fig. |2] shows this is indeed the 
case, by comparing the results of a simulation with zero potential with the analytic 
free-particle solution. 

It will prove convenient to express units in terms of G, h, and a. The basic unit 
of distance is then a, the basic unit of time is (cr^/ (G/i))^/^, and the basic unit of mass 
is {h'^ /{Ga)y^^. These units exhibit the scale-invariance of the Schrodinger-Newton 
equation: a given solution remains a solution if the distance is multiplied by A, time is 
multiplied by A^^^, and mass is multiplied by A~^/^, (and tp is renormalized accordingly). 
For an example of a particular physical scale, if a = 4.47 x 10~^m, as in [1], then time 
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Figure 2. The Li norm of the relative error of a free particle wave function is shown, 
at two different resolutions, as a function of time. The scaling is such that for the 
expected 6th order accuracy the two curves should coincide. 



is given in units of 6.33 seconds and mass is given in units of 3.33 x 10 ^''kg (2 x 10^° 
u). 

Since the computational domain must be finite, the outer boundary is located at 
an arbitrary position but far enough away that remains negligible there throughout 
most of the evolution. The boundary condition is chosen to be Neumann. This choice is 
arbitrary, due to the following machination. As in [20j , in order suppress reflections, the 
potential is modified with the addition of an imaginary tanh term, constructed so as to 
be of order 1 within ~ a of the outer boundary. This term is only non-negligible near 
the outer boundary, where it dampens waves of wavelength less than ~ cr, effectively 
simulating their propagation beyond the computational domain. 

Some additional nontrivial computations are performed for purposes of analysis. 
The expectation value of the Hamiltonian, [H) = {tp\H\ilj), is computed in order to show 
its approach to an eigenenergy as ip approaches a stationary state, with the formula. 



+ m$| 

2m 



dV (12) 



where the differentiation and integration is performed numerically as described above. 
In addition, it is useful to find the peak of r^l-j/'p. This is done by locating, among 
the discrete radii, the point at which r^|'?/'p is maximum, then fitting this and the 
neighboring points to a polynomial. 
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3. Results 

For m > l.U{h^/{Ga)y/^ the wave packet generally fragmented in two. A fraction of 
the probability, with wavelength of order fi^/ (Gm^), appeared to propagate outwards to 
infinity. The rest remained mostly localized in a finite region. Fig. |3] shows the amount 
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Figure 3. Total probability in the computational domain is shown, as a function of 
time, for three different outer boundary locations R. Most of the probability is clearly 
localized within a finite region, and asymptotes to a constant. 

of probability remaining in the computational domain during the course of a simulation, 
for m = 2.09{h'^ / {Ga)y^^ and three different domain sizes. It is apparent that most of 
the probability is highly localized, and constant. 

This remnant wave function appeared to eventually approach a stationary state, 
about which it performed damped oscillations. Fig. H] shows several examples 
of the position of the peak of r^l^/'p oscillating about a central value, with the 
amplitude of oscillation gradually diminishing over time. For l.U{h^ /(Ga))^/^ < m < 
l.53{h^ / {Ga)Y^^ the wave packet initially expanded, then eventually contracted, before 
approaching the apparent stationarity. For m > 1.53{fi^ / {Ga)y^^ most of the wave 
packet contracted immediately. Fig. |5] shows r^l'?/'^ late in the evolution for several 
different masses, after the wave function had become nearly quiescent. The radius of 
the peak of r^l-i/'p was found to scale roughly with /(Gm^), as does the width of the 
MPT groundstate solution found in [12]. The relationship of these wave packets with 
the MPT groundstate solution was confirmed by computing (iplHlip), which was found 
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Figure 4. Radius of tlie peak of r^ji/jp versus time for four different masses. In each 
case the position of the peak appears to perform damped oscillations about a central 
value. 

to approach p^Eq (Fig. |6]), as in [15], where p is the fraction of probabihty that did not 
disperse away and 

Eo = -O.ieSG^mV/i^ (13) 

is the MPT groundstate energy found in [121 [131 [Hi- 
lt is tempting to draw an analogy between the quasistationary final state described 
above and the groundstate of an atom or other conventional potential well. However, 
the p^ factor above represents an important distinction. As first discovered by Harrison 
et al., a time-dependent Schrodinger-Newton solution can evolve towards what might 
be described as a fractional groundstate. Moroz, Penrose, and Tod originally found 
that the energy of the lowest energy stationary state, assuming unit normalization, was 
given by Eq. f lT3|) . However Harrison et al. discovered cases of a wave function that 
partly dissipates away (as also demonstrated in Fig. [3]) while the rest evolves towards 
ip'f^ = p'^ijjQ{pr,t), where p is the remaining probability and ipo represents the MPT 
groundstate with unit normalization. (Note that in general if iplr, t) is a solution of the 
Schrodinger-Newton equation, then so is fi'^iplfir).) The norm of the wave function ip^ 
can be easily verified to be p, and it can similarly be shown {iIjq\H\iIjq) = p^Eq. 

For m < 1.14(^7(^0-))^/^ the entire wave packet appeared to expand indefinitely. 
Fig. [7] shows the trajectory of the peak of r^lT/'p for m = 1.14(^^/(Gcr))^/^ and 
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Figure 5. The probability density shown late in evolution for four different masses, 
all starting from the same initial state, also shown. Each of the "final" states is 
approximately stationary, oscillating only slightly about these amplitudes. The radius 
of the peak scales as fi^ /(Gim?). 



m = 1. 17 {fi^ / {Ga)y^^ . The latter reaches a maximum radius and turns back while 
the former continues expanding. 

This threshold is consistent with computation of the Newtonian escape velocity. 
Noting from Eq. ([3]) that the effective gravitational mass density is proportional to the 
probability density, ^escape = \/ 2 Gnipeak /''"peak , where rripeak is defined as the effective 
gravitational mass within r < rpeak and rpeak is the radius of the peak. Fig. |H] compares 
the velocity of the peak with its escape velocity for four different masses. The velocity 
for m = 1.14(^^/(Gcr))^/^ is just above the escape velocity (after some evolution), while 
the velocity for m = 1. 17 {fi^ / {Ga)y^^ is below the escape velocity, in agreement with 
Fig. [71 The velocities for masses below m = 1.14(/i^/ (Ga)Y^^ were all found to be above 
the escape velocity and their corresponding trajectories appeared to escape to infinity. 
The velocities for masses above m = LUih"^ / {G(r)y/^ were all found to be below the 
escape velocity and their corresponding trajectories were clearly bound. 

The phenomena described above is also roughly consistent with expectations based 
on a "force" -balancing argument of the type given in pLj. The initial acceleration due 
to quantum diffusion for the peak of r^l-^p, where \ip\'^ is Gaussian, and neglecting 
the gravitational potential, is rpeak = ^^/("^^^peak)- Meanwhile the gravitational 
acceleration of the peak is — Gmpeak/^pgak- a Gaussian wave function, mpeak ~ 
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Figure 6. {H) computed during tlie evolution of four masses. Each approaches the 
groundstate energy Eq = — O.lGSG^m^/Ti^, rescaled by the cube of the probability. 

0.43m, so these accelerations are equal and opposite when m ~ 1.32{h'^ / {Ga)Y^^ . Of 
course this argument is only heuristic, but its estimation comes close to the midpoint 
in between the ingoing and outgoing solutions. 

4. Discussion 

In summary, for a wave packet of width a and m > 1.14(^^/(Gcr))^/^ Schrodinger- 
Newton solutions evolve towards the MPT groundstate solution. Since the width of 
the MPT groundstate solution scales as h'^/(Gm^), it can be quite narrow relative 
to the initial wave packet and the interim evolution resembles the traditional picture 
of the "collapse" of the wave function to a Dirac delta function, albeit through 
smooth evolution, and with some fraction of the probability dispersing away. And 
for m < 1.14:{h'^ /{Ga)y^^ the wave function seems to spread to infinity, although more 
slowly than a free particle. 

Schrodinger-Newton theory makes predictions that deviate significantly from 
conventional linear quantum mechanics when the width of a wave packet is large relative 
to /(Gm'^). Unfortunately, the most massive particle used in a quantum diffraction 
experiment to date, a "perfluoroalkylated nanosphere" [21], yields h^/{Gnr') ~ 10^^ 
meters. A more reasonable wave packet width of 1 micron, which is close to what can 
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Figure 7. The radius of tlie peak of r^|?/;|^ vs. t for to = l.lA{h^ / {Ga)Y/^ and 
TO = 1.17(?i^/(G(T))^/'^. The latter is bound, the former appears to propagate outwards 
without bound. 



be achieved in a laboratory, would require a mass of over 10^ u to undergo a Schrodinger- 
Newton collapse. Matter-wave interferometry may not yet be technologically feasible in 
this case, but state-of-the-art experiments are expected to reach 10^ u for silica spheres 
laser-cooled in an optical cavity [22], or even 10^ u for gold clusters beamed through 
diffraction gratings [23]. 
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